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We show that integro-differential generalized Langevin and non-Markovian master equations can 
be transformed into larger sets of ordinary differential equations. On the basis of this transformation 
we develop a numerical method for solving such integro-differential equations. Physically motivated 
example calculations are performed to demonstrate the accuracy and convergence of the method. 



O 
o 

C 

(N 



> : 

00 . 
00 ■ 

o ■ 
vo : 
o . 

m ■ 
O ■ 

^ : 

9 L,: 
^— > ■ 

G ■ 
cd 

3 : 
cr 



INTRODUCTION 

Generalized Langevin equations (GLE) Q and non- 
Markovian master equations 0, 0, 0], which arise in the 
treatment of systems interacting with environmental de- 
grees of freedom, often have an integro-differential form. 
Unlike ordinary differential equations which can be read- 
ily solved using Runge-Kutta, Predictor-Corrector and 
other well known numerical schemes^ there are no gen- 
eral methods for solving equations of integro-differential 
type. Here we show that these integro-differential equa- 
tions can be converted to ordinary-differential equations 
at the expense of introducing a new time variable which 
is treated as if it is of spatial type. [Similar schemes are 
employed to numerically solve the Schrodinger equation 
for time-dependent Hamiltonians0 and as analytical 
tools 7]. There is also some resemblence to schemes for 
solving intego-differential equations of viscoelasticityj^.] 
We then develop a numerical method based on this exact 
transformation and show that it can be used to accurately 
solve a variety of physically motivated examples. 

Neglecting inhomogeneous terms resulting from noise, 
for simplicity, the generalized Langevin equations^ for 
position q(t) and momentum p(t) of a damped oscillator 
in one dimension can be expressed in the form 



dq(t)/dt = p(t)/m 



(1) 



dp(t)/dt = -muj 2 q(t)- / j(t,t')p(t') dt' (2) 



where m and u) are the mass and frequency of the os- 
cillator and j(t, t') is the memory function. Defining a 
space-like time variable u and a function 



X (t,u) =/(«) / 7 (t + u, *>(*') dt', (3) 

J — oo 

it can be verified by direct substitution that p(t) and 
x{t,u) satisfy the following ordinary differential equa- 
tions 



dp(t)/dt 
d\{t, u)/dt 



-muj 2 q(t)- X (t,0) 

/(u)7(* + «,*M*H 
/'(«) 



dx(t,u) 
du 



/(«) 



X(t, u) 



(4) 



(5) 



Here we have introduced a differentiable damping func- 
tion f(u) (with /(0) = 1) which plays a useful role in 
the numerical scheme we will introduce to solve the or- 
dinary differential equations Q, and ©■ [Note that 
/'(«) = df(u)/du.] 

Neglecting inhomogeneous terms, non-Markovian mas- 
ter equations 0,0,0 can be written in the form 



dp(t)/dt = -i[H(t),p(t)] 



K(t,t')p(t') dt' (6) 



where p(t) is the time-evolving reduced density matrix 
of the subsystem, H (t) is an effective Hamiltonian, and 
K(t,t') is a memory operator. [We employ units such 
that h — 1.] Defining an operator 



X(t, «) = /(«)/ K(t + u,t')p(t') dt', 



(7) 



it can be verified by direct substitution that p(t) and 
x(t, u) satisfy ordinary differential equations 



dp(t)/dt 
dx(t, u)/dt 



i[H(t),p(t)]- X (t,0) (8) 
dx{t, u) 



f(u)K(t + u,t)p{t) 
/'(«) 



du 



X(t,u). 



(9) 



Here f(u) is again a differentiable damping function such 
that /(0) = 1. 

Thus, the integro-differential Langevin equations JT}- 
|2"|l can be expressed in the ordinary differential forms Q 
and (@J-© and the integro-differential master equation 
(|o]) can be expressed as the ordinary differential equa- 
tions l|HJl-lj!|Jl. To exploit these transformed equations as 
a practical numerical scheme we must discretize the u 
variable on a grid of points so that the number of ordi- 
nary differential equations is finite. Once this is achieved 
the ordinary differential equations can be solved using 
standard techniques 0|. We use an eighth order Runge- 
Kutta routine|lO| in our calculations. 

To minimize the number of grid points we choose a 
damping function f(u) which decreases rapidly with u. 
In the calculations reported here we used f(u) — e~ gu . 
In practice fewer grid points are needed for positive u 
than for negative u, and we found that the points Uj — 



2 



(— n+l+j)Au for j = 1, . . . , n worked well when we chose 
/ = int(.338n). Here u n — lAu is the largest positive 
u value. While accurate solutions can be obtained for 
almost any non-zero value of g we found the most rapid 
convergence when values were optimized for the type of 
equation. Hence, g is specified differently below for each 
type of equation. To complete the numerical method 
we need a representation of the partial derivative with 
respect to u on the grid. This could be performed via 
fast fourier transform techniques We chose instead to 
employ a matrix representation 



d_ 



3,k 



{j - k)Au 



(10) 



which is known as the sinc-DVR (discrete variable 
representation) 9j. A discrete variable representation 
(DVR) is a complete set of basis functions, associ- 
ated with a specific grid of points, in which functions 
of the variable are diagonal and derivatives have sim- 
ple matrix representations^. DVRs are often used in 
multi-dimensional quantum mechanical scattering theory 
calculations^. In the sinc-DVR[9| , which is associated 
with an equidistantly spaced grid on (— oo,oo), partial 
derivatives can thus be evaluated with a sum 



fdX(t,', 



du 



fc=i 



(j - k)Au 



X(t,u k ) (11) 



for any function or operator X(t,u). In our calculations 
we chose Au to equal the time interval At between output 
from the Runge-Kutta routine. 

We now discuss applications of the above numerical 
method to specific models. For the generalized Langevin 

FIG. 1: Memory functions W(t) plotted against time. 




equation we chose an initial value problem (i.e. j(t, t') = 
for t < t' and j(t, t') = W(t-t') for t > t') where W(t) 
has one of the following forms 



W(t) 
W{t) 



1 e 



9 1 



1 9 



(12) 
(13) 



W(t) 
W(t) 



2e- 2t 
3e" 2t 



e 

2.8e" 



-t/2 



(14) 
(15) 



which arc displayed graphically in Figure 1. The solid 
curve is Ijl2|) , the dashed is (|13|l , the short-dashed is (|14|) 
and the dotted is Ijl5|l . These memory functions were 
chosen to roughly represent the various functional forms 
which can occur physically^] and for ease in obtaining 
exact solutions. The constants appearing in equations 
(@J and 10 are chosen as m = 1 and to 2 = 10. Fig- 
ure 2 shows the functional form of the exact solutions 
q(t) (solid curve) and p(t) (dashed), which evolve from 
initial conditions q(0) = 1 and p(0) = .1, for memory 
function (|12l) over a timescale of 20 units with At = .04. 
Solutions for the other memory functions (and the same 
initial conditions) are similar in appearance. These ex- 
act solutions were obtained by expoiting the fact that 
the above memory functions are sums of exponentials 
(i.e. W(t) = J2jLi a,je~ bjt ) from which it follows that 
one may write 



dp(t)/dt = -muj 2 q(t) - '^2a j e- bit y j {t) (16) 

3=1 

dy 3 {t)/dt = e b ^p(t) (17) 

for j = 1, 2, . . ., and solve these ordinary differential equa- 
tions using standard methods. This approach only works 
for memory functions of this type. Approximate solu- 



FIG. 2: Position (solid curve) and momentum (dashed) of a 
damped oscillator. 




tions were obtained using g = 7/[(n — l)Au] 2 . For nega- 
tive u we set W(u) = 

The negative logarithm of the absolute error in q(t), 



e (t) = - log 10 \q(t) - (^approximate^)!, 



(18) 



is shown in Figure 3 plotted against time for the values 
of n indicated in the inset. [The error in p(t) is simi- 
lar.] As n increases e increases (on average) and hence 
the error decreases. The oscillations in e are caused by 
periodic intersections of the two solutions. In practice 
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it is impossible to visually distinguish the two solutions 
when e > 2. Note that after a short transient the er- 
ror (on average) does not increase. This is probably a 
consequence of the linearity of these equations. Some 
decline in accuracy with time should be expected when 
the Langevin equations are non-linear (e.g. a particle in 
a double- well). 



FIG. 3: e(t) for memory function il l 2t . 
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Figure 4 compares the exact solutions for q(t) (solid 
curve) and p(t) (short-dashed) with those obtained using 
our method for n = 150 (dashed and dotted, respectively) 
over a time of 40 units. No disagreement is visible. Con- 
vergence for memory function (|13fl is similar. 



FIG. 4: Comparison of exact and approximate position and 
momentum of a damped oscillator. 




Memory functions Ijl4|l and Ijl5|l which take negative 
values and have long time tails require many grid points 
for convergence. Figure 5 shows the negative logarithm 
(base ten) of the absolute error in q(t) for this case. While 
many grid points are required, high accuracy solutions 
can clearly be obtained using our method. 

For the master equation we chose an initial value prob- 
lem consisting of a dissipative two- level system represent- 
ing a spin interacting with environmental degrees of free- 



dom. If the spin Hamiltonian is H 



J T a z + j3a x and 



FIG. 5: e(t) for memory function 11411 . 




then the equation for the density matrix p{t) is of the 
formli 



dp(t) 
dt 



/3(r x ,p(t)] 



-C / W(t - t>){a 2 xP (t>) + p{t')a 2 x - 2a x p(t')a x } dt' 
Jo 

(19) 

where the sigmas denote Pauli matrices. Parameters 
were set as ui — 1 = (3 and C = .2. We chose to de- 
fine X (t,u) = f Q W{t - t')p{t') dt' which differs some- 
what from the general definition employed in J2J. The 
transformed equations are then 



dp(t) 
dt 

d X (t,u) 
dt 



Pa x ,p(t)}-2C{ X (t,0) 

a xX (t,0)a x } (20) 
dx(t,u) 



e~ 9U W{u)p(t) 
%gu X (t,u). 



du 



(21) 



Theory predicts that the memory function W(t) for this 
problem is approximately gaussian in form|4|. How- 
ever, we were unable to obtain an exact solution of 
the master equation for this case[ll|. Instead we ap- 
proximate the gaussian via the similar function W{t) — 
14 e -7.4t _ 13e -st_ Exact so i utions f or 

(a z )(t) = Tr{a z p(t)} = Pll (t) - p 00 (t) (solid - curve) (22) 
(a x )(t) = Tr{a x p(t)} = p w (t) + p 01 (t) (dashed) (23) 
((J y )(t) = Tr{a y p(t)} = i(pxo(t) - Poi(t)) (short - dashed) 

(24) 

and initial conditions (<j z }(0) — 1 and (a x )(Q) = = 
(<J y )(0) were obtained in the same way as for the gen- 
eralized Langevin equations and are plotted vs time in 
Figure 6. For the approximate method we used g = 
ll/[(n — I) Au} 2 and for negative u we set W(u) = W(\u\). 
From Figure 7 where we plot 



the coupling to the environment is proportional to a x 



e(t) = -log w \(a z )(t)-(a z ) 



approximate 



(t)\ (25) 
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FIG. 6: Spin x (solid curve), y (dashed) and z (short-dashed) 
components. 
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tions. Because this transformation is exact we expect 
that the method will also work for equations not consid- 
ered in this manuscript. It should be possible to obtain 
accurate solutions for such equations via the following 
steps. First find an approximation of the memory func- 
tion or operator which will allow exact solutions to be 
obtained. Optimize the numerical method by finding the 
best g for the model equations. Finally, apply the nu- 
merical method to the original equations and look for 
convergence of the solutions with increasing n. 

The author gratefully acknowledges the financial sup- 
port of the Natural Sciences and Engineering Research 
Council of Canada. 



against time we see that convergence of the numerical 
method is very rapid for these equations. [Similar accu- 
racies are achieved for (a x ) and (<J y ).] 

FIG. 7: e{t) for (cr z > 



'50' 
75' 
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Thus, we have shown that accurate solutions of 
integro-differential equations can be obtained via trans- 
formation to a larger set of ordinary differential equa- 
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